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A global network of very sensitive large-scale laser interferometric gravitational wave detectors is 
projected to be in operation by around the turn of the century. The network is anticipated to bring 
a range of new astrophysical information — relating to neutron stars, black holes, and the very early 
universe — and also new fundamental physics information, relating to the nature of gravity in the 
strongly nonlinear regime for example. This information is borne by gravitational waves that will 
typically be very much weaker than the level of intrinsic strain noise in the detectors. Sophisticated 
signal extraction methods will therefore be required to analyze the network's data. Here, the noisy 
output of a single laser interferometric detector is examined. A gravitational wave is assumed to 
have been detected in the data. This paper is concerned only with the subsequent problem of 
parameter estimation. Specifically, we investigate theoretical lower bounds on the minimum mean- 
square errors (MSE) associated with measuring the parameters that characterize the waveform. The 
pre-merger inspiral waveform generated by an orbiting system of neutron stars or black holes is ideal 
for this study. Monte Carlo measurements of the parameters of noisy inspiral waveforms have been 
performed elsewhere, and the results must now confront statistical signal processing theory. Three 
theoretical lower bounds on parameter estimation accuracy are considered here: the Cramer-Rao 
bound (CRB); the Weiss- Weinstein bound (WWB); and the Ziv-Zakai bound (ZZB). The CRB is 
the simplest and most well-known of these bounds, but suffers from a number of limitations. It 
has been applied a number of times already to bound gravitational wave measurement errors. The 
WWB and ZZB on the other hand are computationally less simple, and we apply them here to 
gravitational wave parameter estimation for the very first time. The CRB is known as a local 
bound because it assumes that the parameters one seeks to estimate are deterministic, and provides 
bounds on their MSE for every possible set of intrinsic parameter values. The WWB and ZZB 
are known as global (Bayesian) bounds because they assume that the parameter set is random, of 
known prior distribution. They bound the global MSE averaged over this prior distribution. We 
first set up a model problem in order to develop intuition about the conditions under which global 
bounds are more appropriate than their local counterparts. Then we obtain the WWB and ZZB for 
the Newtonian-form of the coalescing binary waveform, and compare them with published analytic 
CRB and numerical Monte Carlo results. At large signal-to-noise ratio (SNR), we find that the 
theoretical bounds are all identical and are attained by the Monte-Carlo results. As SNR gradually 
drops below ~ 10, the WWB and ZZB are both found to provide increasingly tighter lower bounds 
than the CRB. However, at these levels of moderate SNR, there is a significant departure between 
all the bounds and the numerical Monte Carlo results. We argue that the WWB and ZZB are 
probably within a few percent of the theoretical minimum MSE attainable for this problem. The 
implication is that the maximum likelihood method of parameter estimation used by the Monte 
Carlo simulations is not the optimal estimator for this problem at low-to- moderate SNR. In fact, 
it is well-known that the optimal parameter estimator is the conditional mean estimator. This, 
unfortunately, is notoriously difficult to compute in general. We therefore advance a strategy for 
implementing this method efficiently, as a post-processor to the maximum likelihood estimator, in 
order to achieve improved accuracy in parameter estimation. 
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I. INTRODUCTION 



A network of very sensitive instruments is presently being assembled across the globe with a view to directly 
detecting cosmic gravitational radiation on Earth. Although detection is the initial goal of the network, it must 
ultimately function as an astrophysical observatory. In this capacity it will study a range of sources of gravitational 
waves, such as neutron stars, black holes, and the early universe. The network will also provide data through which 
to learn new and fundamental physics about the gravitational field. Arguably, broadband instruments based upon 
the method of laser interferometry offer the best long-term potential for making astrophysical observations. Several 
such instruments have been funded and are being constructed at the present time. They include: LIGO, the US 4 km 
arm-length Laser Interferometric Gravitational Wave Observatory [Q; VIRGO, a French/Italian project to construct 
a 3km arm-length interferometer GEO600, a UK/German effort to build a 600m arm-length interferometer [||; 
TAMA, a Japanese project for the construction of a 300 m arm-length interferometer Joint observations between 
these interferometers are scheduled for soon after the turn of the century. All of the instruments will continue to 
improve their sensitivity incrementally for many years thereafter. 

Interferometers are intrinsically noisy instruments and the presence of noise will mask the identity of all but the 
very strongest incident gravitational waves. This requires a very careful and sophisticated processing of the data, in 
order to extract the valuable information that is borne by the waves |^ . Let us consider the data acquired by only 
a single detector in the network. Suppose that a detection of a gravitational wave of some assumed known form has 
already been made in the instrument's noisy output. Our focus of attention here is the subsequent measurement of 
the one-or-more parameters that characterize the waveform. Specifically, we would like to compute the theoretical 
minimum mean-square errors with which the parameters can in principle be measured. This will clearly impact on 
the astrophysical inferences that can then be drawn as a result of the observation. 

The particular gravitational waveform that we will focus on in this paper is the chirp of gravitational radia- 
tion that precedes the coalescence of a compact binary system comprising of neutron stars (NS) and/or black 
holes (BH). Coalescing binaries are the most promising sources of gravitational waves in the long run for the 
LIGO/VIRGO/GEO600/TAMA detectors (^-0]. As radiation reaction drives the stars through a slow inspiral phase 
just prior to coalescence, the binary generates a very clean gravitational wave signal that is amenable to theoretical 
modeling (see [^,0-jlO| and references therein). LIGO/ VIRGO anticipate observing the last few minutes of NS-NS 
inspiral, during which the gravitational waves oscillate through ~ 16 000 cycles as their frequency sweeps through 
the visibility bands of the detectors. The coalescing binary event rate predictions are subject to gross uncertainties. 
However, it is not unreasonable that there may be a few NS-NS, NS-BH, and BH-BH mergers out to a distance of 
~ 200 Mpc in a period of a year ||ll|-0 . 

Inspiraling binary gravitational waves are encoded with a rich suite of physical and astrophysical information. 
This ra nges from tests of general relativity [^16 1^, through measurements of neutron star and black hole mass and 



spin 1^,^ 18-2^, to new and independent inferences about the value of the cosmological parameters |^l|-^. The 
informational content of binaries has driven gravitational wave theorists to focus much of their collective effort on 
waveform calculations for inspiraling binaries. In tandem with this, there has been considerable study of algorithms 
for analyzing noisy gravitational wave data to extract the waveform information. 

The main goal of this paper is to re-asses the issue of information extraction with respect to observations of coalescing 
binaries by an interferometric gravitational wave detector. Our motivation is the result of a recent confrontation 
between the theoretical Cramer- Rao low bound (CRB) on parameter estimation errors [25[2^ for coalescing binary 
waveforms, with real measurement errors based upon application of the maximum likelihood (ML) method of parameter 
estimation to simulated data sets |2^ . The ML errors were found to depart significantly from the CRB at moderate 
signal-to-noise ratio's (SNR) of around ~ 8. This implies one of the following: (i) the CRB is a weak lower bound at 
this SNR, in which case a tighter theoretical bound would be desirable; (ii) the CRB is a tight bound and the ML 
method is not the optimal estimator, in which case a more refined parameter estimation method would be desirable; 
(iii) the CRB is both weak and the ML method is not optimal, in which case it may be desirable to seek an improved 
theoretical lower bound and an improved estimator. 

In an attempt to discriminate between these options, we apply new theoretical bounds on measurement accuracy 
to the problem. Specifically, we investigate the Weiss- Weinstein bound and the Ziv-Zakai bound These 
are Bayesian bounds that are much more versatile than the more familiar CRB. Although a little more difficult to 
compute than the CRB, they can often be considerably tighter. This has been demonstrated for a range of parameter 
estimation problems in radar and sonar |30| . 

The paper is organized as follows. Sec. II is an overview of parameter estimation, emphasizing the differences 
between local bounds on parameter estimation accuracy of which the CRB is an example, and Bayesian bounds of 
which the Weiss- Weinstein and Ziv-Zakai bounds are examples. The Weiss- Weinstein bound in described in some 
detail in Sec. Ill This is followed in Sec. IV by a detailed description of the Ziv-Zakai bound. In Sec. V, some of the 
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computational issues posed by these bounds are presented. The bounds are apphed, in Sec. VI, to a simple problem 
that illustrates the general superiority of Bayesian bounds over their local counterparts. Then, in Sec. VII, the bounds 
are computed for a Newtonian coalescing binary waveform immersed in Gaussian random noise of spectral density 
characteristic of the first-stage LIGO detectors. We compare the bounds with actual parameter estimation errors that 
were obtained recently from a Monte-Carlo experiment designed to test the maximum-likelihood method. The main 
results are discussed in Sec. VIII and some pointers for future work are given. 



II. SUMMARY OF PARAMETER ESTIMATION 



A common problem in many fields, ranging from radar and sonar through to geophysics and astronomy, is to seek 
estimates for the set of parameters characterizing a waveform that is corrupted by additive Gaussian noise. Consider 
the following such observation 

x{t) = s{t;e) + n{t), \t\<T/2, (2.1) 

and assume the signal s{t]6) is a known function of time for all values of the parameter vector 9, and n{t) is a 
zero-mean Gaussian random noise process. Some measurement algorithm, which we do not need to specify in detail 
yet, is applied to h{t) in order to extract 9. The algorithm produces an estimate 9, with associated error e = 9 9. 
A statistical summary of the performance of the algorithm is contained in the error covariance matrix, given by 
IZ = (ee"*"), where (•) denotes expectation. The diagonal elements of TZ are the mean-square errors (MSE's) on each 
individual parameter, while the off-diagonal elements represent their cross-covariances. In order to benchmark the 
performance of any practical parameter estimator, one would like to know the theoretical minimum TZ for the problem 
at hand. If the signal parameters enter nonlinearly, then a closed-form expression for this cannot be found, and in 
general it is prohibitively difficult to compute numerically. A schematic picture of how the MSE will depend on SNR 
in a nonlinear parameter estimation problem is given in Fig. 1. In the small error or asymptotic region, characterized 
by high SNR, estimation errors are small. In the ambiguity region, where SNR is moderate, large errors occur. When 
SNR is very small, the observations provide little information and the MSE is close to that obtained simply from the 
prior knowledge about the problem. In this paper we will be concerned with bounds that are able to characterize 
performance in the asymptotic and ambiguity regions. These bounds generally fall into one of two classes: local 
bounds or global Bayesian bounds. We will describe their main features in the next Section. 



A. Local bounds 



The formulation of local bounds is based on the premise that the unknown parameters one seeks to measure are 
deterministic quantities. The bounds are local in the sense that they are placed on the MSE's for each different 
possible value of the intrinsic parameter vector. Local bounds have two serious limitations. First, they are restricted 
in application to estimators that are unbiased. In practise, biased estimation is often unavoidable. If the space of a 
parameter is finite, for example, then an unbiased estimator of it does not exist. Secondly, local bounds are unable 
to incorporate any prior information that one might have about the parameters. The Cramer-Rao bound (CRB) is a 
familiar example of a local bound and it therefore has only limited utility. This bound states that for any unbiased 
estimator of a parameter vector 9, based upon noisy observations x, the error covariance matrix must be larger than 
or equal to the inverse of the Fisher information at 9. Thus 

TZ.J ^ m - e,){e, ~ 9,)) > (2.2) 

where is the Fisher information matrix whose elements are given by pq] 



and A(x; 0) is the likelihood ratio 



A(x;.) = «. (2.4) 
p(x|0) 



The CRB is not difficult to compute and it is widely invoked. In particular, it has been used almost exclusively to 
bound the measurement errors on the parameters of gravitational wave signals M,|9|J3|-j33] . Moreover, it can be 
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proved that the CRB is asymptotically attained by the maximum-likelihood (ML) method of parameter estimation 
|p5| . As gravitational wave observations of coalescing binary signals will be rare, and the majority of detections will be 
made at only moderate SNR's, it is unlikely that the asymptotic conditions will be met in practise. Similar situations 
exist in the fields of radar and sonar, and here alternative bounds to the CRB have been considered. The Barankin 
bound, for example, is a local bound that can be much tighter than the CRB [3^ ]. However, it is considerably more 
difficult to compute as it requires maximization over a number of free variables. Also, being a local bound, it still only 
applies to unbiased estimators and is unable to incorporate prior information about the parameters if this information 
is available. 



Rather than treating the unknown parameters as deterministic quantities, Bayesian bounds treat them as random 
variables with known prior distributions. These bounds are global in the sense that they bound MSB's on each of the 
parameters, averaged over their prior distributions. In contrast to their local counterparts, Bayesian bounds are not 
restricted in application to unbiased estimators. In fact, they lower bound the performance of any estimator. Also, 
unlike local bounds, they easily incorporate any prior information about the parameters. It is straightforward to form 
a Bayesian version of the CRB, by simply replacing the conditional probability density p(x|0) with a joint probability 
density p(x, 6) using Bayes theorem. However the Bayesian CRB is subject to a stringent regularity condition: it 
requires the prior probability density function of the parameters, p{6) to be twice differentiable. In the common case 
of parameters that have uniform priors, this regularity condition is obviously not met. 

Another example of a Bayesian bound is the conditional mean estimation bound (CMB) |^6|. In fact this is not really 
a bound, since it can be attained by the conditional mean estimator (CME). This estimator achieves the minimum 
MSE and provides the benchmark against which the performance of other estimators should be compared. The CME 
of a scalar parameter 9, based upon noisy observations x, is given by 



In the context of gravitational wave parameter estimation, the CME has been referred to as the Kallianpur-Striebel or 
nonlinear filter |37| . Unfortunately, the CMB is prohibitively difficult to compute for all but the simplest of problems. 
It generally requires multi-dimensional integrations to be performed numerically over the prior parameter space. 

The complexity of the CMB has motivated the formulation of two further important Bayesian bounds — the Weiss- 
Weinstein bound and the Ziv-Zakai bound. These trade off some of the computational complexity of the CMB, and 
yet are only apparently a little less tight. In a range of applications the bounds have demonstrated this utility 
We now describe each of these bounds in some detail in the following two Sections. 



The Bayesian form of the CRB and the CMB discussed in the previous Section belong to a general class of Bayesian 
bounds that Weiss and Weinstein were able to derive from the Schwarz inequality [^8| . An outline of their derivation 
is given here. The WWB is a member of this general class, but it is free from the problems that limit the CRB and 
CMB. It therefore is of much more general utility than the latter two bounds. A version of the WWB for the case 
of a single parameter is obtained below. A statement of the multiple-parameter generalization of the scalar bound 
follows. 



A lower bound on the error in estimating a scalar parameter 9, based upon noisy observations x, is sought. Let 
p(x, 9) denote the joint probability density of x and 9. Weiss and Weinstein introduced a function ■(/'(x, 9) such that 



B. Bayesian bounds 




(2.5) 



III. WEISS- WEINSTEIN BOUND 



A. Single parameter 




Vx. 



(3.1) 



Since, for any real-valued measurable function ,g(x). 
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/oo pOC 
dxg(x) / d0V(x,0Mx,0) =0, (3.2) 
-oo J —oo 

the condition in ( |3.l| ) implies that ■0(^5 ^) is orthogonal to any transformation of the data x. Subtracting (6'?/'(x, 6)) 
from both sides above and then applying Schwarz's inequality to the left side, results in the following 

^- (^2(x,g)) ■ ^^-^^ 

As this inequality is valid for any ^(x), it sets a lower bound on the mean-square error in estimating 6 from observation 
of X. To un derline this point we will replace (/(x) in the following by ^(x). Note that the lower bound set by the right 
side of (3_^) is independent of the estimator (i.e. is absent of ff(x)). Recall that local bounds, such as the CRB, do 

not share this property of Bayesian bounds: they apply only to estimators that are unbiased. 

It is simple to see that the Bayesian form of the CRB and the CMB are special cases of the general bound in (^.3|). 
Consider choosing the function '(/'(x, 6) as follows: 

(3.4) 



This choice satisfies the orthogonality condition (3.1) and generates the Bayesian CRB. Similarly, the selection 

V'(x, 61) = - (^Ix) (3.5) 

leads to the CMB. 

In their quest for a less restrictive bound than the CRB and the CMB, Weiss and Weinstein were led to consider a 
different choice for ip(x, 6). They proposed the following 

7/'(x, 61) = L''(x; 9 + d,e)~ L^-''{k; 9-6,6), (3.6) 

where r and S are arbitrary real-valued scalars and L(x; 9i, 62) is the likelihood ratio 

Li.;eM-'-^y (3.7) 
P(x, 6*2) 



This choi ce fo r V'(x, 9) satisfies the orthogonality condition ( |3.l| ) for all combinations of 6 and < r < 1. Substitution 
into Eq. (3.3) generates the WWB on the mean-square error, e^, in the estimation of 9: 



g2 > (5^ exp[2?7(r,(5)] 

- exp[ry(2r, S)] + exp['q{2 - 2r, -S)] - 2 exp[?7(r, 26)] ' ^ ' ^ 



where 

r]{r,5) = In {L''{x;9 + 6,9)) 

= In p'{9 + 5)p^-'{9) I j p'{Ti\9 + 5)pi-'^(x|6l)dx| d6' 

= ln/ /(6' + %i-'^(6l)exp[^(r;6'-f^(5,6l)]d6l. (3.9) 

Several comments are pertinent here. First, note that the integration with respect to 9 is performed over the region 
Q = {9 : p{9) > 0} in order to avoid singularities. Second, the bound reduces to the Bayesian version of the CRB for 
6 0. Third, the term /i(r; 9 + 6,9) is a familiar one in information and communications theory: it is known as the 
semi-invariant moment generating function and used to bound the probability of error in binary hypothesis testing 
problems |26[ . We shall meet it again in our discussion of the Ziv-Zakai bound in the next Section. 

In order to remove one of the degrees of freedom, the WWB is usually computed for r = 1/2. It then reduces to 

, ^^exp[27y(l/2,^)] 
- 2{l-exp[ry(l/2,25)]}- ^ ' ' 

The variable 6 that enters the bound is usually referred to as a test point. The optimal value for 6 is the one that 
generates the maximum bound. This value may be other than 6^0, for which the WWB reduces to the CRB as we 
remarked earlier. 
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Weiss and Weinstein have shown how to generahze the single test point bound (3.1C) to incorporate multiple test 
points. Consider a vector of N test points S = {6i, 62, ■ ■ ■ , 5n). The corresponding multiple-test point WWB is 



> uQ-^u' , (3.11) 

where the elements of the vector u are 

= S^ , (3.12) 

and the elements of the matrix Q are 



,exp[77(l/2, - 6j)] - exp[?7(l/2, 6, + 6^)] 
exp[277(l/2,50] 



(3.13) 



In order to evaluate Eq. (3.11), a matrix of dimension equal to the number of test points has to be inverted numerically. 
This imposes a practical restriction on the number of test points. However, as we shall see later, the WWB fortunately 
appears to converge quickly with increasing N. 



B. Multiple parameters 

Since the coalescing binary waveform is characterized by more than one parameter, we shall require a multiple 
parameter version of the WWB. Consider a vector of M parameters, 9 = {6i, . . . ,0m)- The WWB on the error 
covariance matrix TZ is obtained in a similar fashion to the single parameter bound. The result is, 

n>ng-^n'^. (3.14) 

The elements of the matrix H are the M x N test points in the multi-dimensional parameter space. The N x N 
matrix Q has elements given by 

C _ . exp[r;(l/2, S, - S,)] - exp[77(l/2, S, + 6,)] 
exp[7y(l/2,^,)] cxp[77(l/2,5,)] ' 

where di is the i'th test point in the parameter space and 



r;(l/2,J,(5j) =ln ^ ^p(6/ + S,)p{e + S^) p{x\9 + S^)p{x\9 + J,)dx| d9 . 



(3.16) 



Again, several comments are pertinent. First, integration with respect to 9 is over the region O ~ {9 : p{9) > 0}. 
Second, for a non-singular bound there must be at least M linearly independent test points. Finally, ( |3.14 ) reduces 



to the Bayesian CRB upon setting Ti. — 6J, where J is the identity matrix and the scalar 5 — > 0. We will turn to the 
practical issues involved in the computation of the scalar and vector parameter versions of the WWB in Sec. V. 



IV. ZIV-ZAKAI BOUND 

The theoretical foundation of the Ziv-Zakai bound (ZZB) is somewhat different to the WWB there do not 

appear to be any formal theoretical links. In common with the WWB however, the ZZB places a fundamental lower 
bound on the performance of any parameter estimator. We present a simple derivation here, for the case of a single 
parameter with a uniform prior distribution. The multiple-parameter extension of the bound for arbitrary priors is 
also presented. 



A. Single parameter 

As a concrete example, suppose that an estimate of the difference between the arrival time of a gravitational wave 
at two separated detectors is required. Let us denote this parameter by 9. Now ask what is the probability of 
making a correct decision between two possible values, </> and (/) -I- A, of this parameter. The likelihood ratio test 
(LRT) is the optimal decision scheme that produces the minimum probability of error. Instead of the LRT, consider a 
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simpler suboptimal decision scheme in which a decision is made in favour of the "nearest-neighbour" to some arbitrary 
estimate, 9, oiO. Thus, 



Decide Ho-.e^cj) if e<(j)+ — , 



Decide Hi 



A if 



> 



A 
"2" 



(4.1) 



If the two hypothesized delays are equally likely to occur, which is physically most reasonable, then the suboptimal 
decision scheme has a probability of error given by 

A 



. + A) = iP^>0+- 



A 



(4.2) 



Clearly if P„ 



• + A) is the minimum probability of error, associated with the LRT, then 



A) < ( ^ > ^ 



(4.3) 



where e = 9 — 9 denotes the estimation error. Now, suppose that 9 is uniformly distributed on [—T,T]. In t his s pecific 
example, T would represents the gravitational wave travel time between the two detectors. The inequality (4.3) holds 
good for any (j) and A, in particular combinations of </> and A such that + A e [— T, T], or 



-T<6<T-A, 



< A < 2T . 



Integrating Eq. (E^) with respect to (j) over [— T, T — A] gives 



T-A 



+ A)d0 < 



P i\^\> 



A 



which can equivalently be expressed as 



T-A 



A)d(t> <T H 



where 



iJ(A)^^ J%{\e\>A\ 



(4.4) 



(4.5) 



(4.6) 



(4.7) 



Note that (4.6) is only useful for A < 2T, since f or A > 2T the integral is negative and therefore zero is a better 
bound. The next step is to multiply both sides of (4.6) by A/T and integrate with respect to A over [0, 2T]. Noting 
that 



2T 



A2d{iJ(A)} 



(4.8) 



is the mean-square error in the estimation of 9 when the latter has a uniform prior distribution in [— T, T], the 
integration yields: 



2r 



2T 



AdA 



T-A 



Pn 



A)dA, 



(4.9) 



which is the ZZB in its simplest incarnation. Bellini and Tartara ]39| have remarked that i?(A) is a non-increasing 
function of A and suggested that the bound might be tightened by applying a 'valley-filling' function to the left-side 
of (4.6). Denoting this function by V[-], the Bellini- Tartara version of the ZZB is 
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^ - 2T 



1 i'2T j-T-A 



dA. (4.10) 



The bound generalizes in a straightforward manner for an arbitrary prior, to give 

e'>/ ^-vlf [p(0)+p(0 + A)]-^'min('^,0 + A)d4dA. (4.11) 



Jo 2 l^J-oo J 

Although there is generally no closed form expression for Pminl'/', + A)j tight lower bounds exist p6[ . 

B. Multiple parameters 

The ZZB has only recently been extended to vector random parameters with arbitrary prior distributions 
Consider an Af -dimensional vector random variable, 6, with prior pdf p{6). As before, let 6 be an estimate of 9 
produced by any estimator, e the estimation error, and TZ = (ee"'") the error covariance matrix. Then the following 
lower bound on sl^TZel for any M-dimcnsional vector a has been obtained: 

f°° A 

a"^7^a> / — -y'dA, (4.12) 
Jo 2 

where the valley filling fimction V' is now defined as 

T/' = l/|max J{p{<f>)+ p{<j) + S))-P„U<t>,(P + S)d4>'j , (4.13) 

with the maximum referred to S. The bound is generated, as for the single parameter case, via an inequality between 
the probability of error in a suboptimal decision rule and the minimum probability of error associated with an LRT. 
However, one has now to decide between one of two possible values cf) or (p + S for the parameter vector under 
investigation. The suboptimal decision rule is then: 

T- T A 

Decide Hq : 6 = 4> if a'0>a'0+y, 
Decide Hi : 6 ^ 4> + 6 if sJe<sJ(f>+^. (4.14) 



The hyperplane 



aT0 = aV+y, (4.15) 



separating the two decision regions, passes through the midpoint of the line connecting (p and cf)+6 and is perpendicular 
to the a axis. A decision is made in favour of the hypothesis that is on the same side of the separating hyperplane as 



the estimate 9. The tightest bound in ( [4. 12]) is achieved by maximization over the vector S, subject to the constraint 
a"''<5 = A. This constraint does not determine the vector d uniquely. In order to satisfy it, S must be composed of a 
fixed component along the a axis, and an arbitrary component orthogonal to a. That is 

<5 = ^^a + b, (4.16) 

II a ll-' 

where 

a"^b = 0, (4.17) 

and there are M —1 degrees of freedom in choosing 8 via the vector b. Simply setting b = results in hypotheses that 
are separated by the smallest Euclidean distance. However, this does not necessarily guarantee the largest probability 
of error. A maximization over 5 can improve the bound. 
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V. COMPUTATIONAL ISSUES 



In this Section the CRB, WWB, and ZZB are reduced to forms that are appropriate for calculating error bounds 
on the parameters of a signal that is immersed in a background of Gaussian random noise. 

The signal waveform s(t, 9) is parameterized by a vector, 6, for which an estimate is sought. Noisy measurements 
of the signal are obtained as follows: 

x{t)=.s{t,e)+n{t), -|<i<|, (5.1) 
where n{t) is the noise, assumed Gaussian with a known spectral density denoted by Sn{f)- 



A. The Cramer-Rao Bound 



The CRB was defined in (2.2) to be the inverse of the Fisher information at 6. The latter is a matrix of second 
derivatives of the likelihood ratio of an observation with respect to 6. In the case of stationary Gaussian noise, the 
matrix elements reduce to 



7-1 — 



ds 



(5.2) 



where (si|s2) denotes the inner product between two signals, si{t) and S2(t). In terms of the signal's Fourier trans- 
forms, si(/) and S2(/), and the spectral density of the noise, Sn{f), the inner product can be expressed as 

(.rM = 2y^ SM) ^'-'^ 

The integral is a measure of the degr ee o f 'overlap' between the two signals, and in radar applications it is often 
termed the ambiguity function. Thus, ( ^.2| ) can be interpreted as the local curvature of the signal ambiguity function 
around its maximum. Numerical integration is generally required to compute the elements of the Fisher information 
matrix. It is then straightforward to perform the inversion and obtain the CRB. The diagonal elements of the inverted 
Fisher matrix are the Cramer-Rao bounds on the variances of each of the signal's parameters. 



B. Weiss- Weinstein Bound 



The c alculation of the WWB relies upon evaluating the semi-invariant moment generating function ri[l/2;6i,Sj) 
in ( 3.16 ). It is not difficult to show (see Appendix in |2^ for details) that, for stationary and Gaussian noise, this 
function can be reduced to 



where 



77(1/2; S,,Sj) = In C{S„ S,) + /i(l/2; - Sj) , 



(5.4) 



(5.5) 



and the region of integration is = {0 : p{0) > 0}. In Eq. (5.4), the first term embodies the prior information about 
the parameters. Consider the single test-point version of the WWB for a single parameter having a uniform prior on 
the interval [—D/2, D/2]. The integral is simple to evaluate and it yields 



C(<5)=ln (1-1^ 



(5.6) 



After some algebra, the second term in (5^) reduces to 



m(1/2;^. 



(5.7) 
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where = (s|s) is the squared ampHtude (energy) signal-to- noise ratio, and 



{s[e + 8, 8,] I s[e]) 

-l{5, - Sj) = ' (5.8) 

{s[0] I s[0]} 

is the normahzed signal ambiguity function. It is often the case in practise that the latter function is independent 
of 6, and then it's calculation is greatly simplified. This will be the case for the examples that are presented later. 
However, numerical integrations are still generally required to compute the signal ambiguity function. Moreover, 
these integrations have to be performed for every set of test point locations in the parameter space. As the WWB 
also requires inversion of a matrix having dimension equal to the number of test points, it is clearly desirable to keep 
the number of test points down to a minimum. An indicator of the number of test points that are required for a 
given problem, and their optimal locations, is the shape of the signal ambiguity function. As we shall see later, for 
the coalescing binary waveform this function has a very well-defined shape. 



C. The Ziv-Zakai Bound 



The main term on which the evaluation of the ZZB in ( 4.12 ) hinges is Pmin, the minimum probability of error in a 
binary detection problem. An exact expression for Pmin exists for the decision problem of discriminating between two 
equally likely signal vectors, si and S2, in a background of Gaussian noise of covariance K. The minimum probability 
of error is then given simply by 



(5.9) 



where d is the normalized distance between the signals 



(s2-Si)TK-l(S2-Si) 



(5.10) 



and 



/2n 



(5.11) 



If the inner products under the square root sign in (5.10) are evaluated, one finds that 



^-^i{l/2;8,-8,), 



(5.12) 



where /i is given by Eq. (5.7). Therefore the calculation of the ZZB, like the WWB, is crucially dependent on the 
shape of the signal ambiguity function. In the case of the WWB this dictates the number of test points and their 
locations in order to achieve a tight bound. In terms of the ZZB, the shape of the signal ambiguity function defines 



a path of integration in (4.12), subject to the constraint under which the integration is evaluated. 



VI. AN ILLUSTRATIVE EXAMPLE 



In this Section we apply the CRB and one of the Bayesian bounds (WWB) to a simple scalar parameter estimation 
problem. Our intention is to investigate the conditions under which Bayesian bounds arc tighter than local bounds. 

A common feature of all the bounds that we have presented is that they depend on the shape of the signal ambiguity 
function, 7, rather than the signal shape, when the background of noise is Gaussian. Often the shape of the signal 
ambiguity function is the same for all underlying values of the signal's parameters, greatly simplifying the calculation 
of the bounds. 

The CRB only probes the shape of the signal ambiguity function around it's maximum. Structure in the ambiguity 
function away from the maximum could be enhance by noise and masquerade as a false peak. This would confound 
a maximum likelihood parameter estimator, and may lead to numerical parameter estimation errors that depart 
significantly from the theoretical GRB. 

While the GRB is 'blind' to the presence of sidebands in the signal ambiguity function, the WWB and ZZB are able 
to capture this structure. In the case of the WWB this is achieved through the test points. As well as probing around 
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the main lobe of the ambiguity function, test points may also be placed around the secondary maxima. Similarly, the 
ZZB is generally tighter than the CRB if a path of integration is selected to traverse all of the predominant lobes in 
7- 

The difference between the CRB and the WWB is best illustrated through an example. We consider two signals, 
si and S2, characterized by a scalar parameter, 9, that enters the signals nonlinearly. The signal waveforms need 
not concern us here, only their ambiguity functions 7,(0, where i — 1,2. We have chosen signals whose 7's are 
independent of and depend only on displacements in 0, i.e. AO. These ambiguity functions are displayed in Figj^. 
The signal si was designed so that 71 has only a single broad maximum. The other signal, S2, has 72 comprising of 
a number of significant secondary lobes. Note also that 71 is actually the "envelope" of 72 . 

We assume that 9 has a uniform prior. However, the region of support of the prior is set much larger than the 
anticipated parameter estimation errors, so that the prior does not actually impact upon the estimation accuracy for 
this problem. 

For si, the WWB was computed by placing N test points uniformly along the lobe of the signal ambiguity function. 
It was found that the resulting bound was not very sensitive to where the test points were placed along the lobe for 
this signal. A variable number of test points (up to 20) were used to study the convergence of the bound. This was 
attained for only 4 test points. For S2, the test points were placed around the main lobe of the ambiguity function 
and also around the principal secondary maxima. In fact only the first three secondary lobes needed to be covered: 
again the WWB exhibited rapid convergence. 

The CRB was also obtained for si and S2 over an identical range of SNR. This was obtained by inverting the Fisher 
matrix ( |5.2| ) for this specific problem. However, the CRB can also be computed in terms of the WWB formalism for 
a single test point, 5, allowing S ^ 0. It should be remarked that in the multiple test-point formulation of the WWB, 
one test point is always forced to have (5^0. This ensures that the WWB reduces to the CRB at large SNR. 

The bounds, e^'^ (i = 1, 2), that we calculated are displayed in the top panel of Fig. ^. For si, the WWB and the 
CRB are in good agreement down to SNR ~ 11. At smaller values of SNR, the WWB is a tighter bound than the 
CRB but not significantly so. This result is not too surprising: 71 has only a single broad maximum and the CRB 
probes the curvature of this lobe. Similarly the WWB probes the structure in the ambiguity lobe, although it is able 
to probe further away from the maximum with respect to the CRB. As the lobe is broad, the WWB is generally a 
little tighter than the CRB. The results for S2 are significantly different. Here, the WWB departs from the CRB at 
a much higher value of SNR, around 20. This is because the CRB is blind to the secondary maxima in 72 that the 
WWB is able to capture through a judicious choice of test points. The discrepancy between the CRB and the WWB 
for this example is striking: a factor of '--^ 5 at SNR =10 and more than an order of magnitude at SNR = 5. The 
WWB falls significantly as SNR increases, while the CRB remains fairly constant, due to the sharp maximum in 72 . 

It is also interesting to compare the behaviour of e^^'' and e^^\ At high SNR the accuracy in the determination of 

9 for S2 is better than for si. This is clearly due to the sharp maximum in 72. At low SNR, we expect e"g^ ^ ejp 
because the ambiguity functions have roughly the same "global" profile. This intuition is borne out in the results of 
the Bayesian analysis (WWB), but not for the local analysis (CRB) where < at all SNR. 

VII. APPLICATION TO COALESCING BINARY PARAMETER ESTIMATION 

We are now in a position to investigate whether Bayesian bounds provide a tighter constraint than the CRB on 
the error covariance matrix for the parameters of a gravitational wave signal generated during the inspiral phase 
of a compact object binary. In the following we use units where G — c — 1. This implies a conversion factor: 
1M„ = 4.926 X 10-S. 



A. signal and noise model 

The coalescing binary inspiral waveform can be cast in the following generic form 

h{t)^A[nf{t)f/^cos[4>{t)]. (7.1) 

where f{t) is the instantaneous gravitational wave frequency, 0(i) the instantaneous phase, and A the amplitude. 
We will consider the phasing of the wave, as well as the amplitude evolution, only up to Newtonian order. In this 
approximation the factor A is a constant, complicated, function of the binary's distance, location in the sky, chirp 
mass M. = vr?-J^'°m'2^ l(ra\ + m2)^l^ (where m\ and m2 are the masses of the compact objects) and the detector's 
antenna pattern |^. It's precise functional form need not concern us further. Frequency and phase read 
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fit) = fa 



8/3 



(7.2) 



and 







-5/3| 




fa 





(7.3) 



where the constant r, sometimes referred to as the chirp time, can be cast in terms of the chirp mass of the binary as 



256 ^ ' 



(7.4) 



The constants fa and $a are respectively the frequency and phase of the signal at the arbitrary time t = ta- The 
waveform is characterized in terms of 3 parameters: ta^ <I>a, and r. The amplitude parameter A ent ers the waveform 
linearly and it will be incorporated later into our definition of signal-to-noise ratio (see Eq. (7.12)). Of course, the 
parameterization of the signal is not unique and one can express h{t) as a function, for example, of the time to 
coalescence and the phase of the wave at this time. In fact, the latter parameterization may have some advantages. 
However, our main goal here is to assess the relative difference between the CRB and Bayesian bounds and so this 
detail of the parameterization is not crucial: we require only consistency in the choice of parameterization of the signal 
in order to compare the bounds. In particular the set of parameters that we are assuming here may not correspond 
to physical ones, and this would be the case if we extended this setting to post- Newtonian waveforms [|^. However, 
there is one crucial feature of the wave parameterization that we adopt here: it produces a signal ambiguity function 
that is independent of the intrinsic values of the parameters and depends only upon their displacements. This fact 
is more transparent if we examine the signal's Fourier transform ft.(/), which in the stationary phase approximation 
reads: 



M/)=AA/-^/6 exp[vI-(/)] , 



(7.5) 



where 



2r 



N = ( ^ ) ft" 



1/2 



(7.6) 



is a normalization constant, and 



*(/) = z^^,(/)A^-zJ; 



represents the parameter vector 



(7.7) 



y={ta,<^a.T) 



(7.1 



and 



^-2 =-1, 
V-a = 2Tif - 



, fa 



-5/3 



(7.9) 



Notice that the signal's paramet ers enter linearly into the phase ( [7.7j) of its Fourier transform. The normalized 
ambiguity function for the signal (7.5) is given by 



7(AA^ 



(7.10) 



where 
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(7.11) 



and the integral is defined over the frequency interval, within the instrument's sensitivity band, spanned by the signal. 
The optimal signal-to-noise ratio reads: 



EE {h\h) = 



(7.12) 



and, therefore, incorporates the amplitude parameter A via the definition of A/", cfr. Eq. (7.6) 



Our model for the noisy spectral density, 5'„(/), is intended to be representative of the performance of the first stage 
LIGO detectors. An analytic fitting formula for this has been presented in |3l], and we utilize it here. Accordingly to 
pTf , we have considered the observational window confined to the frequency interval 40 — 750 Hz. We suppose that 
the final frequency of inspiral is outside the considered bandwidth, so that the integral involved in the definition of 
the inner product (5.3) is evaluated on the same frequency range. 



B. Calculating the bounds 



In order to compare the bounds on parameter estimation errors given by local and global approaches, we computed 
the CRB, WWB and ZZB. The computational steps are described here, focusing particularly on the WWB and the 
ZZB. Our discussion is centered upon the evaluation of bounds for the time-of-arrival parameter, ta (= Ai). Similar 
results apply to <I>a and r. 

The CRB involves the computation of the diagonal elements of the inverse of the Fisher information matrix ( ^.2| ), 
that is: 



£2 = 7-1 • 
'-'it ' 

for the signal (^^) and the parameters ( |7.§| ), we have 



dh 



and, therefore, 



Snif) 



df. 



(7.13) 



(7.14) 



(7.15) 



The evaluation of (|7.13| ) is straightforward from Eqs. (7^) and ( 7.15 ) and has been thoroughly studied in many 



papers ||l4|lJ3|J3a-BJ where further details can be found 



The WWB involves the computation of Eq. ( |3.14 ) and therefore of the M x N matrix H and the N x N matrix Q, 
where M is the number of parameters (3 in this problem) and N is the number of test points {Sj; j=l, . . . , N}. The 
test points are now given explicitly by AXj (the lower and upper indices labelling the test points and the parameters, 
respectively). The elements of the matrix Ti are therefore Hi.j = AA^. 

We will assume here that X'^ has a uniform prior distribution with a region of support that is much larger than the 
anticipa ted e rrors on th e parameters. Therefore the prior will not impact on the calculation of the WWB. 
Eqs. (3.15) and (5.7) read now 



Gjk 



and 



^exp[p2^(AA- - AA-)/4] - exp[p2^(AA- + AA-)/4] 
exp[p2^(AAp/4] exp[p2-y(AAp/4] 



77(1/2; AA^, AAD = Ml/2; AA^ - AAD . 



(7.16) 



(7.17) 



A crucial issue for a reliable computation of the bound is the placing of test points: (i) in order to get the CRB in the 
limit of high SNR, the elements Hi.j for j = 1,2,3 (notice that j runs from 1 to N) have been chosen accordingly to 
Wj/j = (5yctAA", with AA" <C 1 (here S^a is the Kronecker symbol); indeed, we always probe the primary peak of the 
ambiguity function, as does the CRB; (ii) in order to get the tightest possible bound at low SNR, we placed the other 
test points (Tiuj for j > 3) along the maxima of the ambiguity function, using the test-case problem presented in Sec. 
VI as a guidehne. In Fig. || we show 7(AA''), in the plane (A^a, At), maximized with respect to A$a (we already 
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know that the regions where the ambiguity function is small do not contribute significantly to the result). The plot 
is enlightening as 7 consists on a long, sharp ridge and clearly indicates that the test points need to be spread along 
that curve. We placed up to 25 points (about the maximum permitted by the numerical routines implemented for 
the matrix inversion) with different choices of their separation and distance from the origin. We noted, in fact, that 
with only 4 almost equally spaced points (spacing ~ 16 ms), the result did not change significantly, in agreement with 

what found in the toy problem. 

The evaluati on o f the ZZB involves the computation of the integral (4.12), using the minimum probability error 
given by Eqs. and ([sl^ ). As we have stressed before, the strategy of computation replaces here the spreading 

of test points with the selection of the integration path. Our discussion of the evaluation of the WWB indicates that 
the integration has to be performed along the ridge of the ambiguity function shown in Fig. ^, in order to produce 
the tightest bound. Of course, no "valley filling" function was needed, as 7 is a smooth curve free from oscillations 
for the signal that we were studying. We carried out the integration up to a maximum displacement from the origin 
~ 0.2 sec (of the same order of the position of the last test point during the investigation of the WWB), after which 
no appreciable improvement was found. As in the case of the computation of the WWB, the prior probability on Ai 
did not impact on the bound, because the region of support was chosen to be much larger than the anticipated error. 



C. Results 



The bounds calculated following the three different theoretical approaches (CRB, WWB and ZZB) were computed 
for values of p in the relevant range 7 < p < 25. In the same SNR interval we compared these results with those 
obtained by means of numerical simulations that implement maximum likelihood estimators. The root-mean 
square error bounds on the time-of-arrival parameter ta are displayed as a function of p in Fig. ^ All the bounds 
converge to the CRB at p ^ 15, and at this SNR the CRB is attained by the maximum likelihood estimator. At 
smaller values of SNR, the Bayesian bounds deviate from the CRB, providing a slightly tighter result (~ 6% at p = 10 
and ~ 25% at p = 7). This not too severe discrepancy can be explained by the structure of 7: both the local and 
global approaches probe the ambiguity function around its origin, but the Bayesian bound is able to follow it further 
away from the origin (see the discussion in Sec. VI). The behavior of WWB and ZZB were found to be very similar, 
although not exactly equal: the SNR threshold at which they depart from the CRB is p ~ 14, for the ZZB, and 
p ~ 12, for the WWB, but the latter provides a better constraint at low SNRs. This agrees with results from other 
applications of the bound to time-of-arrival parameter estimation problems in radar and in sonar applications (see 1 3^] 
and reference therein). The striking feature of the comparison is given by the maximum likelihood errors obtained in 
numerical experiments: while matching the behaviour of (local and global) lower bounds at high SNR (p ^ 15), they 
produce errors that are dramatically higher than expected theoretically at low SNRs: about 65% at p = 10 and more 
than a factor of 2 at p = 7. 



VIII. CONCLUSIONS 



We have reassessed the issue of information extraction with respect to observations of coalescing binaries by inter- 
ferometric gravitational wave detectors. After discussing the main properties of global and local bounds on parameter 
estimation, we applied a set of these bounds to a gravitational wave parameter estimation problem. In particular 
we have introduced the Weiss- Weinstein and the Ziv-Zakai bounds. These provide fundamental lower limits on the 
mean-square error on the parameters that describe a signal, independently of the actual estimator that is adopted 
in the data-analysis process. In addition these bounds easily incorporate any a priori information that is available 
about the problem and do not suffer from limitations that affect local bounds and the Bayesian version of the Cramer- 
Rao bound. In short, these global bounds can be used to benchmark the performance of any practical information 
extraction technique. 

We have applied the bounds to the case of laser interferometric measurements of waveforms that are characteristic 
of those emitted by inspiraling compact binaries in a background of noise that is characteristic of the performance of 
first-stage detectors. Comparisons between the Cramer-Rao, Weiss- Weinstein and Ziv-Zakai bounds on the MSE and 
actual maximum likelihood errors obtained by numerical experiments, over a wide range of SNR, show that: (i) at 
high signal-to-noise ratios (SNR ^ 15) all the approaches converge to the same value of the MSE. In this regime one 
can regard the Fisher information matrix as a simple and reliable tool to compute "realistic" bounds on estimation 
errors. Maximum likelihood methods are probably adequate to extract astrophysical information from noisy data 
at these SNR's (although a definitive statement is premature, pending a detailed analysis applied to more general 
waveforms); (ii) at low signal-to- noise ratios (SNR ^ 10, in which most of the events are likely to be recorded during 
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the first years of operation of the detectors) the WWB and ZZB produce a more stringent constraint (~ 25% at 
/9 = 7) on the MSE with respect to the CRB, indicating that the latter can underestimate the errors in this regime. 
Perhaps more seriously, all the bounds are about two times smaller than the errors that are obtained in the numerical 
experiments. 

This analysis suggests that maximum likelihood techniques need to be refined, or complemented, in order to attain 
the lowest possible value of the errors. Our study of toy-problems and the results from previous investigations in 
different fields suggests that these are within a few percent of those predicted theoretically via the WWB and ZZB. A 
first attempt toward the understanding of the outcome of numerical experiments and the performances of maximum 
likelihood estimators has been recently reported in [Q. We are currently exploring the possibility of implementing 
conditional-mean estimators to improve parameter estimation accuracy, and will report on this in a forthcoming 
paper p3| . The conditional mean-estimator is generally intractable to implement as it requires the computation of a 
multi-dimensional integral over the space of all the parameters (in realistic cases more than ten), even though suitable 
strategies to reduce the amount of computation have been proposed and successfully tested in (simple) cases Q. 
Hierarchical strategies combining maximum likelihood and non-linear filtering could also speed up the process. 

Finally, it is important to underscore the point that the form in which we have presented the WWB and ZZB is 
completely general. It can be applied to parameter estimation problems for other kinds of signals {e.g. pulsars) and 
other instruments and/or arrays of detectors. 
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FIG. 1. Schematic representation of the behaviour of the mean square error as a function of the signal-to-noise ratio in a 
typical non-linear estimation problem. 



FIG. 2. This diagram illustrates the link between the structure of the ambiguity function and the accuracy of parameter 
estimation. The lower panel displays the ambiguity function for two different one-parameter signals, as a function of the 
displacement from the true parameter value: si (bold line) was chosen to have an ambiguity function with one broad 
maximum; S2 (dashed line) was designed to have a more structured ambiguity function with many secondary side-lobes. For 
the special signals considered here, both ambiguity functions depend only on AS, and not on the actual value 9. In the upper 
panel we display the results of applying the Cramer-Rao and the Weiss- Weinstein theory to bound the mean-square error ee in 
the estimation of the signal's parameter as a function of the signal-to-noise ratio. See text for further details. 



FIG. 3. The ambiguity function, maximized with respect to the phase-of-arrival $„, for the signal (7.5) as a function of the 
time-of-arrival and of the chirp time. 

FIG. 4. Comparison between local and global theoretical bounds on the time-of-arrival error with actual parameter estima- 
tion errors obtained by applying the maximum likelihood method to simulated data (bold line: CRB; unfilled circles: maximum 
of WWB and ZZB; full circles: maximum likelihood error). 
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